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ABSTRACT 



A method an apparatus uses a variational approach for space 
dependent gamut mapping. The method presents a new 
measure for the problem. A solution to the formulation 
according to the method is unique if the gamut of the target 
device is known. A quadratic programming efficient numeri- 
cal solution provides real-time results. 
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SPACE VARYING GAMUT MAPPING 

TECHNICAL FIELD 

[0001] The technical field is color image processing using 
gamut correction. 

BACKGROUND 

[0002] Gamut mapping is a method used to modify a 
representation of a color image to fit into a constrained color 
space of a given rendering medium. A laser-jet color printer 
that attempts to reproduce a color image on regular paper 
would have to map the photographed picture colors in a 
given color range, also known as the image "color gamut," 
into the given printer/page color gamut. Conventional gamut 
mapping methods involve a pixel by pixel mapping (usually 
a pre-defined look-up table) and ignore the spatial color 
configuration. More recently, spatial dependent approaches 
were proposed for gamut mapping. However, these solutions 
are either based on heurestic assumptions or involve a high 
computational cost. An example of such a method is dis- 
cussed in "Color Gamut Mapping Based on a Perceptual 
Image Difference Measure/' S. Nakauchi, et al., Color 
Research and Application, Vol. 24, pp. 280-291 (1999). 

[0003] Another method relies on preserving the magnitude 
of gradients in an original image, while projecting the 
gradients onto the target gamut as a constraint. This multi- 
scale property is achieved by sampling the image around 
each pixel with exponentially increasing sampling intervals 
while the sampling is done along vertical and horizontal 
axes. The method involves modulating an measure for 
image difference by human contrast sensitivity functions. 
The method uses a model in which the contrast sensitivity 
function is a linear combination of three spatial band-pass 
filters H ls H 2 , H 3l given in the spatial-frequency domain (or 
hi, b 2 , h 3 , as their corresponding spatial filters), as shown in 
FIG. 1. 

[0004] For gamut mapping of the image Uq in the CIELAB 
(which stands for Commission International de FEclairage 
(CIE) with L, A, and B (more commonly L, a, b) respec- 
tively representing converted signals for cyan (C), magenta 
(M) and yellow (Y)) space, the method minimizes the 
functional 



f=t ce\La.b\ Jn 



[0005] subject to {u\ u e , u b }e9. 

[0006] In Equation (1), is the filter corresponding to 
the spectral channel c e{L» a, b} the i e{l, 2, 3} 'contrast 
sensitivity' mod, Q is the image domain, and 8 is the target 
gamut. Note that a total of nine filters H a c are involved, three 
for each spectral channel and a total of three spectral 
channels. 

[0007] The filters are modeled by shifted Gaussians. 
H lC is not shifted, and thus, is also Gaussian, while h 2 c 
and h 3 c are Gaussians modulated by two sine functions with 
different periods. A graphical analysis of h 2 c and h 3 ° as 
shown in FIG. 2 reveals that h£ and h 3 ° approximate the 
derivative operator at different scales. These two gradient 



approximation operators are denoted by V o1 ° and Voi c , 
Note that any band pass filter can be considered as a version 
of a derivative operator. Furthermore, one possible extension 
of the ID derivative to 2D is the gradient. Thus, the 
minimization of the Equation (1) functional is similar to 
minimizing the following functional for each channel sepa- 
rately. 

f l^*(«-«o)t 2 + |v^( U -«o)l 2 + |v^{«-« 0 )l 2 rfn (2) 



[0008] Roughly speaking, the first term corresponds to the 
S -CIELAB perceptual measure, while the next two terms 
capture the need for matching the image variations at two 
selected scales that were determined by human perception 
models. One technical difficulty of the spatial filters corre- 
sponding to Equation (1) is their large numerical support. 

[0009] Another simple spatial-spectral measure for human 
color perception was proposed in "A Spatial Extension of 
CIELAB for Digital Color Image Reproduction," Zhang et 
al., Proceedings of the SID Symposium, Vol. 27, pp. 731-34, 
(1996), The 'S-CIELAB' defines a spatial-spectral measure 
for human color perception by a composition of spatial 
band-pass linear filters in the opponent color space followed 
by the CIELAB Euclidean perceptual color. 

SUMMARY 

[0010] Gamut mapping is a method to modify represen- 
tation of a color image to fit into a constrained color space 
of a given rendering medium. A laser-jet color printer that 
attempts to reproduce a color image on regular paper would 
have to map the photographed picture colors in a given color 
range, also known as the image "color gamut," into the given 
printer/page color gamut. 

[0011] A method an apparatus uses a variational approach 
for space dependent gamut mapping. The method presents a 
new measure for the problem. A solution to the formulation 
according to the method is unique if the gamut of the target 
device is known. A quadratic programming efficient numeri- 
cal solution provides real-time results. 

[0012] In an embodiment, the method comprises receiving 
an input image, converting color representations of an image 
pixel set to produce a corresponding electrical values set, 
applying a space varying algorithm to the electrical values 
set to produce a color-mapped value set, and reconverting 
the color-mapped value set to an output image. The space 
varying algorithm minimizes a variational problem repre- 
sented by 

[0013] subject to u e6, where Q is a support of the 
image, a is a non-negative real number, D=g*(u-u 0 ), 
and g is a normalized Gaussian kernel with zero 
mean and a small variance a. 

[0014] The method may be implemented in a variety of 
image capture reproduction devices, including cameras and 
printers. The method may be embodied on a computer- 
readable medium, such as a CD-ROM, for example. 
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DESCRIPTION OF THE DRAWINGS 

[0015] The detailed description will refer to the following 
drawings, in which like numerals refer to like objects, and in 
which: 

[0016] FIG. 1 is a qualitative description of filters mod- 
eling the human contrast sensitivity functions in the spectral 
domain; 

[0017] FIG. 2 illustrates a shifted Gaussian; 

[0018] FIG. 3 is a block diagram of an apparatus that uses 
a variational approach for gamut mapping; 

[0019] FIGS. 4-7 illustrate algorithms and processes used 
with the apparatus of FIG. 3; 

[0020] FIG. 8 illustrates an example of a behavior of an 
algorithm used by the apparatus of FIG. 3; 

[0021] FIG. 9 is a flowchart of an algorithm illustrating an 
alternative operation of the apparatus of FIG. 3; and 

[0022] FIG. 10 is a block diagram of a computer system 
that implements the algorithm of FIG. 4. 

DETAILED DESCRIPTION 

[0023] Gamut mapping is a method used to modify a 
representation of a color image to fit into a constrained color 
space of a given rendering medium. A laser-jet color printer 
that attempts to reproduce a color image on regular paper 
would have to map the photographed picture colors in a 
given color range, also known as the image "color gamut," 
into the given printer/page color gamut, Conventional gamut 
mapping methods involve a pixel by pixel mapping (usually 
a pre-defined look-up table) and ignore the spatial color 
configuration. More recently, spatial dependent approaches 
were proposed for gamut mapping. However, these solutions 
are either based on heurestic assumptions or involve a high 
computational cost. 

[0024] The gamut mapping problem is related to the 
Retinex problem of illumination compensation and dynamic 
range compression. The basic Retinex problem is: How to 
estimate the reflectance image from the given acquired 
image? A reasonable optical model of the acquired image S 
asserts that it is a multiplication of the reflectance R and the 
illumination L images. Where the reflectance image is a 
hypothetic image that would have been measured if every 
visible surface had been illuminated by a unit valued white 
illumination source, and the illumination image is the actual 
illumination shaded on surfaces in the scene. In the log 
domain: 

s-n-l 

[0025] where s, r, and 1 are the respective logarithms of S, 
R, and L. Since surface patches can not reflect more fight 
than has been shaded on them, R<l-*r<0. Thus, in an ideal 
situation, image r<0, which is perceptually similar to s. For 
the Retinex, an additional physically motivated constraint is 
provided, namely, that the illumination image 1-s-r is 
smooth, i.e. the gradient |Vl|-|V(r-s)| is small. But this is 
just another way to say that the features of r are similar to 
those of s, since the illumination is assumed not to create 
perceptual features in s. In the gamut mapping problem an 
image u 0 exists, and the problem is to estimate an image u 
e8, which is not only perceptually similar to u 0 , but also has 
similar perceptual features as u 0 . 



[0026] A good measure of image deviation captures the 
perceptual difference between the initial, Uq and the final, u, 
images. This is modeled by 

D-g m (u-u£ (3) 
[0027] where g may be a normalized Gaussian kernel with 
zero mean and a small variance o. This model is good for 
small deviations. However, for large deviations it should be 
elaborated to account for possible perceptual feature differ- 
ences, which may be modeled by the difference of gradients, 
which due to linearity, turns out to be the gradient of 
Equation (3) 

VZ>=Vk*(H-«o)>g*(V«-Vuo) (4) 

[0028] The proposed measure yields the functional 

E(u)~ f{D 2 + a\VD\ 2 )dil ^ 
Jci 

[0029] subject to u e0, 

[0030] which should be minimized. 

[0031] Taking a first variation of Equation (5) with respect 
to u yields the Euler-Lagrange equation: 

^ = g * (adiv{Vg * (u - u 0 )J -g * (« - uo)) = 0, (6) 



[0032] Reformulating Equation (6) as a gradient descent 
flow for u, provides the following minimization scheme: 

[du (7) 
— = ag*&D -g*D 

< at 

. sj. ued 

[0033] The functional (Equation (5)) and the resulting 
minimization scheme are both Euclidean invariant in the 
image plane. They are thus both translation and rotation 
invariant. As the parameter a goes to zero, the S-C1ELAB 
model is approximated, while for effective a the result is a 
proper extension to the perceptual measures, with an effi- 
cient numerical implementation. 

[0034] To provide a numerical solution, recall that an 
artificial time parameter t was added to the image u(x, y), 
which now reads u(x, y, t). A first step is to discretize the 
Eueler-Lagrange gradient descent equation, by first taking a 
simple forward explicit approximation for the t derivative, 



=ag*&D-g*D 



[0035] where x=dt and u D (x,y)«*u(x,y; nx). 

[0036] Next, the space derivatives are solved. Let 
Uij^-uCihjhjnt), where uniform spatial spacings are assumed 
in the x and y directions of size h. Using central derivatives 
in space, 
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[0037] and the same for the y direction. Using the relation 
g*D xx (g*u)=g 3C *g x *u, the algorithm next computes the ker- 
nels D 2 =g 3t *g x +g y *g y =D x g*D x g+D y g*D y g. The explicit 
approximation reads: 

[0038] subject to the constraint u^" €9. 

[0039] To speed up convergence, a standard coarse to fine 
pyramidal approach may be preferred. For example, an 
original image S may be composed of a set of 512x512 
pixels. The image S may then be averaged down (sub- 
sampled or decimated) by taking, for example, every 4 
adjacent pixels and replacing the block of four pixels with a 
single pixel having a value that represents an average of the 
four pixels. This results in an image set that is 256x256 
pixels. The sub-sampling can be continued, yielding a sub- 
sequent reduced pixel set, or layer, comprising 128x128 
pixels. This process may continue, with each resulting layer 
being an increasingly coarser version of the original image 
S. The final layer will be a set of only 2x2 pixels, and likely 
will be blurred. The process of sub-sampling or decimating 
to provide increasingly coarser details results in an image 
pyramid. The original image S represent the base of the 
image pyramid at 512x512 pixels, followed by layers of 
256x256 pixels, 128x128 pixels, 64x64 pixels, and so on. 
Using this approach of an image pyramid allows for rapid 
convergence of the solution to the functional of Equation 
(5). The process starts by using the coarsest layer of the 
image pyramid, and running a number of iterations until a 
solution is achieved. The solution is then used to initialize 
the next finer layer of the image pyramid, resulting in less 
computation in order to arrive at a solution. The process 
continues until the finest layer (e.g., the original, full reso- 
lution, image set for S) is reached. As an alternative to the 
image pyramid, a full multi-grid method may be used. 

[0040] The functional of Equation (5) has a Quadratic 
Programming (QP) form, since the penalty term is quadratic 
and the constraint is linear. If the set 0 is convex, the overall 
problem is convex if and only if the Hessian of the func- 
tional Equation (5) is positive definite. In such a case, there 
is a unique local minimum that is also the global solution to 
the problem. In the above embodiment, the Hessian is given 
by g*(l-ctV)*g, which is positive definite for all a>0. Thus, 
for a convex target gamut 9, there exists a unique solution. 

[0041] The above-described functional (Equation (5)) may 
be used in imaging devices to map a gamut of an image to 
the gamut of a device that attempts to reproduce the image. 
FIG. 3 is a block diagram of an apparatus 100 that may be 
used for gamut mapping according to the functional of 
Equation (5). An image capture device 101 receives an input 
original image and produces an output image S 102, which 
is an electrical signal representing a calorimetric value 
image. For example, the capture device 101 may convert at 



least three calorimetric values of each pixel of the original 
image into corresponding electrical signals. The electrical 
signals may indicate the L, a, b values, for example. Other 
calorimetric values may be the XYZ tristimilus value and 
the L, U, V or device dependent RGB values, A gamut 
mapper 105 produces an in-gamut image S106. Finally, a 
rendering module 107 provides a rendered image S"108. 
The rendering module 107 may be implemented as a color 
laser printer. The thus-generated image S"108 may represent 
a best-fit image, given gamut limitations of the device in 
which the image reproduction is to occur. 

[0042] The image S 102 may represent the image as 
sensed by the capture module 103. The gamut mapper 105 
applies an algorithm to extract and map the values of the 
image S 102 into the gamut of the image reproduction device 
107. In particular, the gamut mapper 105 may apply an 
algorithm that solves the problem represented by Equation 
(5), thereby solving the gamut mapping problem and opti- 
mizing the output image S"108. 

[0043] FIG. 4 is a block diagram showing routines of an 
algorithm 120 that may be used to complete the gamut 
mapping function. In an embodiment, the algorithm 120 
provides a numerical solution to the functional of Equation 
(5). 

[0044] The algorithm 120 begins with input block 121, 
where an input to the algorithm 120 is an image S of size [N, 
M], and the two parameters a and p. 

[0045] In initialization block 123, a Gaussian pyramid of 
the image s is computed. The thus-constructed pyramid 
contains p resolution layers (from fine (1) to coarse (p)) with 
the current resolution layer (k) set to p. In block 125, T k 
iterations of a gradient descent algorithm are applied to the 
resolution layer, until all resolution layers are checked, 
block 127. In block 129, the next resolution layer is updated. 
When all resolution layers are processed, the result is the 
final output of the algorithm 120. 

[0046] FIG. 5 is a block diagram showing the initializa- 
tion routine 123 in more detail. In block 131, a counter is 
initialized with M. In block 132, the Gaussian pyramid is 
constructed by smoothing the image with a kernel, such as 
the kernel k PYR : 



KPYR = 
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[0047] In block 133, the image is decimated by, for 
example, a 2: 1 ratio. The process is repeated (block 134) and 
the counter is incremented (block 135) p times where p<lg 2 
(min (M, N)). This process produces a sequence of images 
{S] £ } k . 1 p conventionally named the "Gaussian Pyramid of 
S." The image S A is the original image S, and S p is an image 
with the coarsest resolution for the Gaussian pyramid. 

[0048] In block 137, the process sets k=p, i.e., starts at the 
coarsest resolution layer k, and sets the initial condition 
L^-max {S p }. 
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[0049] FIG. 6 is a block diagram of an embodiment of the 
processing main routine 125. The routine 125 starts at block 
141, where the new resolution layer is initialized and T k and 
a k are set (e.g., T k =K*T 0 ). 

[0050] Then, for j=l, ...,T k . 

[0051] In block 143, the routine 125 calculates the gradi- 
ent: 

[0052] where Ax is the convolution of each of the color 
planes of x with Kl^j,: 



0 1 0 

1 -2 1 
0 1 0 



and or* is, for example a k = a *2 2( *" I) 



[0053] In block 145, the routine 125 calculates |^sd : 



Vnsd - 



[0054] In block 147, the routine 125 completes the gradi- 
ent descent iteration: 

[0055] Where n 0 is a constant; for example ^ o -0.8. 

[0056] In block 148, the result is projected onto the 
constraints: 

[0057] Where Proj e (x) is a projection of x into the gamut 
6. 

[0058] Id block 149, if j*T K , processing returns to block 
143. Otherwise, the routine 125 ends. 

[0059] Returning to FIG. 4, in block 127, if k>l, the result 
is up scaled (2:1 ratio) by pixel replication into the new 
Lq, the initialization for the next resolution k-1 layer. The 
resolution layer is updated k=k-l, and the algorithm pro- 
ceeds by going again to block 125. If k«l, the result L^, is 
the final output of the algorithm 120. 

[0060] FIG. 7 is a block diagram of an alternative embodi- 
ment of an algorithm routine, denoted as 120a. The routine 
120a may be applied to solve the equation: 



[0061] subject to the constraint u y D €0. 

[0062] The routine 120a begins with block 151, where the 
kernels D 2 and g, and the counter k are initialized. 

[0063] In block 153, the routine 120a calculates: 
[0064] In block 155, the routine 120a calculates: 



[0065] In block 157, the routine 120a performs the gra- 
dient steps to solve u^ 1 as noted above. 

[0066] In block 158, the routine 120a determines if k>l, 
and if so, processing moves to block 159, where k is set to 
k-1. Processing then returns to block 153. Otherwise, the 
routine 120a ends. 

[0067] The penalty function as shown in Equation (5) 
tends to create halos in the resulting image. FIG. 8 explains 
the origin of those halos through a one dimensional example. 
FIG. 8 shows an original signal 160 that is outside of gamut 
162 (delimited in the gray range by dotted lines 163 and 
164). Projecting the signal 160 onto the gamut 162 will 
result in a constant value, denoted as "High" (163) and loss 
of all detail. FIG. 8 also shows the result of processing the 
original signal 160. Dashed line 165 represents the result of 
scaling the signal 160 into the allowed range of the gamut 
162. All the details are preserved, but with a smaller con- 
trast. As opposed to point operations, the space dependent 
approach represented by Equation (5) yields a signal 166 
(the solid line) that preserves the details with high contrast. 
However, halos occur near strong edges 167, which means 
that near the edges 167 there is a slow transition from low 
to high values. 

[0068] In order to avoid these phenomena, the penalty 
term may be modified, and robust estimation tools may be 
used. The original penalty term in Equation (5) may be 
replaced by: 



(6) 



[0069] which for p A (x)=p 2 (x)«x 2 coincides with Equation 
(5). If the function p(x) grows slower than x 2 as x—»oo, 
behavior near strong edges impr oves. Good candidates for 

p(x) are p(x)«|x| or p(x)=^l+x 2 . 

[0070] A different and simpler (linear) approach with 
similar robust behavior involves solving the original Equa- 
tion (5) twice, with two different values of a. A solution with 
a small a may be denoted as u 8mftn and a solution that 
corresponds to the high value of a may be denoted as u high . 
The solution u amall has smaller contrast at areas with small 
details, yet has almost no halos. On the other hand, the 
solution u high preserves the small details, but at the expense 
of strong halo effects. By averaging these two results (u 8mal , 
and u^J in a spatially adaptive way provides a simple, 
improved solution. The improved solution is therefore: 

[0071] The weight w[kj] should be close to one near 
strong edges, and close to zero in relatively smooth regions. 
In an embodiment, 



1 

1 + $V**«ol 2 



[0072] provides a robust estimation. 

[0073] Halo problems have been recently dealt with in 
relation to Dynamic Range Compression. Solutions pro- 
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posed included anisotropic diffusion and robust filtering. 
The halo-related solutions described herein are solutions to 
the same halo problem. 

[0074] FIG. 9 is a flowchart showing an operation 180 of 
the gamut mapper 105. The process begins in start block 
185. In converter block 195, the raw image S is converted 
from electrical signals having arbitrary values (calorimetric 
values) to locally denned color descriptors. 

[0075] In block 200, the converted image signal is deci- 
mated to form an image pyramid. A 2:1 decimation scheme 
may be used, for example. Subsequent steps in the operation 
180 begin with the coarsest resolution layer of the image 
pyramid using a space varying algorithm. In block 205, the 
gamut mapper 105, using a space varying algorithm such as 
that represented by Equation (5) calculates the image devia- 
tion for the coarsest resolution layer. This process may 
involve calculating a first variation (block 210), determining 
a gradient descent flow (block 215), and solving the result- 
ing gradient subject to a constraint (block 220). In block 225, 
a determination is made if the current resolution layer is the 
last (finest) resolution layer. If additional layers remain for 
processing, the process 180 returns to block 205. Otherwise, 
the process 180 moves to block 230, where local color 
descriptors are converted to the required output format. In 
block 245 the process 180 ends. 

[0076] The above-described space-varying algorithm 120, 
as represented by Equation (5), for example may be 
executed using a general purpose computer system 250 as 
shown in FIG. 10. The computer system 250 includes an 
input device 260, a computer unit 270, a display 280, a 
keyboard 290, and a storage device 300. The storage device 
300 may be a CD-ROM drive, for example. Program code 
instructions for executing the algorithms 120, 120a and 180 
may be stored on a CD-ROM 305, or other computer 
readable memory, which is read by the storage device 300. 
The program code read out of the CD-ROM 305 are 
executed by a CPU of the computer unit 270. The CPU may 
perform the processing shown by the flowchart of FIG. 9. 
The algorithms 120, 120a, and 180 also may be performed 
in a camera or printer hardware. 

In the claims: 

1. A method for gamut mapping of an input image using 
a space varying algorithm, comprising: 

receiving the input image; 

converting color representations of an image pixel set to 
produce a corresponding electrical values set; 

applying the space varying algorithm to the electrical 
values set to produce a color-mapped value set; and 

reconverting the color-mapped value set to an output 
image. 

2. The method of claim 1, wherein the space varying 
algorithm minimizes a variational problem represented by: 



subject to u e8, wherein Q is a support of the input image, 
a is a non-negative real number, D-g^u-uJ, g is a nor- 



malized Gaussian kernel with zero mean and a small vari- 
ance a, u 0 is the input image, and u is the output image. 

3. The method of claim 2, further comprising: 

solving the variational problem at a high value of a; 
solving the variational problem at a low value of a; and 
averaging the solutions. 

4. The method of claim 3, wherein the step of averaging 
the solutions comprises using a spatially adaptive weighting 
scheme, comprising: 

wherein the weight w[kj], comprises: 

w[k, j] = - — — I -r , and 

wherein p is a non-negative real number. 

5. The method of claim 2, wherein the variational problem 
is solved according to: 



— =ag*&D-g*D, subject to u e &. 

6. The method of claim 2, wherein the space varying 
algorithm is solved according to: 

ug wl -V+T(aV-D?J f 
subject to u y D &, wherein 

x-dt, 

Eh,=g*g*(u n -Uo) 
L D =D 2 *(u D -u 0 ) and 

7. The method of claim 1, wherein the space varying 
algorithm minimizes a variational problem represented by: 

E(u)= f(pi(D) + ff ft (|VD|))*fa 

Jo. 

subject to u e8, wherein p 1 and p 2 are scalar functions. 

8. The method of claim 2, further comprising: 

decimating the input image to create one or more reso- 
lution layers, wherein the one or more resolution layers 
comprises an image pyramid; and 

solving the variational problem for each of the one or 
more resolution layers. 

9. The method of claim 1, wherein the method is executed 
in a camera. 

10. The method of claim 1, wherein the method is 
executed in a printer. 

11. A method for color gamut mapping, comprising: 

converting first calorimetric values of an input image to 
second calorimetric values of an output device, wherein 
output values are constrained within a gamut of the 
output device; and 
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using a space varying algorithm that solves an image 
difference problem. 

12. A computer- readable memory for color gamut map- 
ping, comprising an instruction set for executing color 
gamut mapping steps, the steps, comprising: 

converting first calori metric values of an original image to 
second colorimetric values, wherein output values are 
constrained within a gamut of the output device; using 
a space varying algorithm that solves an image differ- 
ence problem; and 

optimizing a solution to the image difference problem. 

13. The computer-readable memory of claim 12, wherein 
the image difference problem is represented by: 



E{u) 



= IV 



+ a\VD\ 2 )dn 



subject to u €0, wherein Q is a support of an input image, a 
is a non-negative real number, D«g*(u-u Q ), g is a normal- 
ized Gaussian kernel with zero mean and small variance or, 
u a is the input image, and u is an output image. 

14. The computer-readable memory of claim 12, wherein 
the instruction set further comprises steps for: 

solving the image difference problem at a high value of a; 

solving the image difference problem at a low value of a; 
and 

averaging the solutions. 

15. The computer-readable memory of claim 14, wherein 
averaging the solutions comprises using a spatially adaptive 
weighting scheme, comprising: 

wherein the weight w[k, j], comprises: 



1 



1 +$Vg*«o|2' 



wherein (5 is a non-negative real number. 
16. The computer-readable memory of claim 12, wherein 
the image difference problem is represented by: 



E(u)= fipi(D)^ a p 2 i\VD\))dCl 
Jci 



wherein p a and p 2 are scalar functions. 

17. The computer-readable memory of claim 12, wherein 
the instruction set further comprises steps for: 

decimating the input image to create one or more reso- 
lution layers, wherein the one or more resolution layers 
comprise an image pyramid; and 

solving the image difference problem for each of the one 
or more resolution layers. 

18. The computer-readable memory of claim 17, wherein 
the instruction set further comprises steps for: 

(a) initializing a first resolution layer; 



(b) calculating a gradient G for the resolution layer, the 
gradient G comprising: 

G=A(u-tO+a K (u-tO, 

wherein Ax is a convolution of each color plane of x with 



0 l 0 

1 -2 1 
0 1 0 



and Oj-Oo^ 20 ^; 

(c) calculating a normalized steepest descent value Lj^Lj. 
i-/4 0 */*nsd*(j» wherein fi 0 is a constant; 

(d) projecting the value onto constraints Proj e (Lj), 
wherein Proj e (x)is a projection of x into a gamut 8; and 

(e) for a subsequent resolution layer, repeating steps 

(b)-(d). 

19. A method for image enhancement using gamut map- 
ping, comprising: 

receiving a input image; 

from the input image, constructing an image pyramid 
having a plurality of resolution layers; 

processing each resolution layer, wherein the processing 
includes completing a gradient iteration, by: 

calculating a gradient G; 

completing a gradient descent iteration; and 

projecting the completed gradient descent iteration 
onto constraints; and 

computing an output image using the processed resolution 
layers. 

20. The method of claim 19, wherein the gradient G, 
comprises: 

G»A (u-u^^ki^ («-«o), 

wherein u is the output image, u Q is the input image, and 
a is a non-negative real number. 

21. The method of claim 19, wherein completing the 
gradient descent iteration comprises calculating: 



Lj = Lj-i -h, Hnsd G, 



and 



wherein j*wsd & a normalized steepest descent parameter,^ 
is a constant, k is a number of resolution layers in the image 
pyramid, and j is a specific resolution layer. 

22. The method of claim 19, wherein projecting the 
completed gradient descent iteration onto the constraints is 
given by: 

wherein Proj e (x) is a projection of x into a gamut 8. 
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23. The method of claim 19, wherein constructing the 
image pyramid, comprises: 

smoothing the input image with a Gaussian kernel; 

decimating the input image; and 

setting initial conductive L Q «max{Sp}, wherein Sp is an 
image with the coarsest resolution layer for the image 
pyramid. 

24. The method of claim 23, wherein the Gaussian kernel, 
comprises: 



Kpyft = 



J_ 1 J_ 
16 8 16 

ill 
8 4 8 

I 1 1 

16 8 16 



E(u)= f(D 2 +a\VD\ 2 )dn, 
Jfi 



subject to u e8, wherein Q is a support of the image, and 
D=g*(u-u Q ), wherein g is a normalized Gaussian ker- 
nel with zero mean and small variance a, u G is the input 
image, u is the output image, and wherein a is a 
non-negative real number. 
26. The method of claim 19, wherein processing each 
resolution layer comprises applying a space varying algo- 
rithm to minimize a variational problem represented by: 



m= f (/>,(z>)+<*ft(ivDi))rfa 



25. The method of claim 19, wherein processing each 
resolution layer further comprises applying a space varying 
algorithm to minimize a variational problem represented by: 



subject to u eG, wherein p x and p 2 are scalar functions. 
27. The method of claim 26, wherein p t and p 2 are c hosen 

from the group comprising p(x)=|x| and p(x)=^l+x 2 . 
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